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Abstract 

A variety of problems in device and materials design require the rapid 
forward modeling of Maxwell's equations in complex micro-structured ma- 
terials. By combining high-order accurate integral equation methods with 
classical multiple scattering theory, we have created an effective simula- 
tion tool for materials consisting of an isotropic background in which are 
dispersed a large number of micro- or nano-scale metallic or dielectric in- 
clusions. 

Keywords: Maxwell equations, multiple scattering, meta-materials, fast 
multipole method 

1 Introduction 

We describe in this paper a simulation method for Maxwell's equations suitable 
for microstructured materials consisting of separated inclusions which are em- 
bedded in a homogeneous background (Fig. [I]). In practice, it is often the case 
that the shape and permittivity of the inclusions are fixed and that one seeks 
to optimize their placement to create a specific electromagnetic response. Each 
new configuration, however, requires the solution of the full Maxwell equations. 
If there are thousands of inclusions in an electrically large region (many wave- 
lengths in size), the calculation is generally too expensive to carry out within a 
design loop. 

In oder to accelerate such calculations, we have coupled complex geometry 
Maxwell solvers with multiple scattering theory. Using the hybrid solver, cal- 
culations such as the one depicted in Fig. [T] require only a few minutes on a 
single CPU, despite the fact that there are a million degrees of freedom needed 
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to describe the full geometry (and there would be orders of magnitude more 
points needed in a finite difference or finite element discretization). 

Our method, which we refer to as fast multi-particle scattering (FMPS), is 
based on a two step procedure. First, we enclose a representative scatterer, 
such as a single pair of gold nanorods, in a sphere S. We then build the scat- 
tering matrix for this nano-structure (described below) using integral equation 
techniques. The solution to the full Maxwell equations can then be obtained 
in geometries with N inclusions (N = 200 in Fig. [I]), by solving the multiple- 
scattering problem where the inclusions have been replaced by their scattering 
matrices. Not only does this reduce the number of degrees of freedom required, 
but we have effectively precomputed the solution operator for each inclusion in 
isolation, so that the linear system we solve by iteration on the multi-sphere 
system is well-conditioned. Further, the fast multipole method (FMM) reduces 
the cost of each iteration from 0{N 2 ) to 0(N log N) and is particularly efficient 
when applied to this problem. 

The principal limitations of the method are (1) that some modest separation 
distance between inclusions is required and (2) that some of the efficiency is 
based on the fact that only a few distinct nanoparticle types are allowed. In 
many experimental settings, both conditions are satisfied. We will return to a 
discussion of these limitations in our concluding remarks. 

2 Maxwell's equations and the Debye-Lorenz- 
Mie formalism 

Working in the frequency domain and assuming a time dependence of e~ lLJt , 
Maxwell's equations in a linear, isotropic material take the form 

VxH w = -iweE tot , (I) 
V x E* ot = iLUfiH to \ 

where E to * and H to * are the total electric and magnetic fields, e is the per- 
mittivity of the medium and /i its permeability. We are mainly interested in 
dielectric inclusions embedded in a background medium, but will consider per- 
fect conductors briefly at the end of this section. The total fields (E*°*,H tot ) 
can be written as the sums of the incident fields (E m , H m ), defined only in the 
exterior region, and scattered fields (E, H) defined in both the inclusions and 
the exterior: 

jjtot = H' m + H. (2) 

It is well-known jTTJ that at dielectric interfaces, the Maxwell equations ([!]) 
are uniquely solvable when supplemented by the the continuity conditions: 

[n x E tot ] =0 =>■ [n x E] = - [n x E m ] 

[n x H*°*] =0 [n x H] = - [n x H m ] (3) 
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Figure 1: Two hundred gold ellipsoid pairs are randomly oriented in the region 
[0,100] x [0,100] x [0,20] and illuminated from above by a plane wave in TE 
polarization. The transmitted z-component of the Poynting vector is plotted on 
planes at z = —4 and z — —8. The wavelength is 2tt so that the particles are 
approximately one wavelength in size, and the region is about 15 x 15x3 wavelengths 
is size. 
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and the Silver-Miiller radiation conditions on the scattered field. The expression 
[n x F] is used to denote the jump in the tangential components of the vector 
field F at a point on the interface. 

2.1 Debye Potentials 

About a century ago, Debye, Lorenz, and Mie [SJ [TSl HO] independently solved 
the problem of scattering from a single sphere by using separation of variables. 
Without entering into the derivation, it is straightforward to verify that 

E(x) = V x V x (xu(x) + iujeV x (xu(x)) 

H(x) = V x V x (xtt(x)) - iuj^V x (x«(x)) (4) 

represent an electromagnetic field, where x denotes the position vector with 
respect to the sphere center, so long as the Debye potentials u, v satisfy the 
scalar Hclmholtz equation 

Au + k 2 u = 0, Av + k 2 v = , 

with Helmholtz parameter (wave number) k 2 = uj 2 t\i. In the exterior of a sphere, 
the Debye potentials u, v can be represented by the multipole expansions 

oo n 

u ext (r,9A) = E E b n , m h n (kr)Y™(9,4>) 

n—0 m— — n 
oo n 

v ext (r,9^) = E E a n ,mh n (kr)Y^(0^) (5) 

n— m— — n 

where (r, 9, </>) are the spherical coordinates of the point x with respect to the 
sphere center, h n (r) is the spherical Hankel function of order n, and Y™(9,<f)) 
is the usual spherical harmonic of order n and degree m. The resulting electro- 
magnetic field then also satisfies the appropriate radiation conditions at infinity. 
In the interior of a sphere, u and v can be represented by the local expansions 

oo n 

u mt (r,9,<f>) = E E d n , m j n {kr)Y™{9,4>) 

n—0 m— — n 
oo n 

v mt (r,9,<f>) = E E Cn,mj»V°r)Y™(p,<l>) (6) 

n—0 m— — n 

where j n (x) is the spherical Bessel function of order n. 
Remark 2.1. To improve readability, we will abbreviate 

oo n 

EE - e 

n— m— — n n,m 
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and the truncated sum 

p n p 

EE as E 



n— m=— n 



If is straightforward to verify that the total number of terms in the truncated 
summation is (p + l) 2 . 

2.2 Single sphere scattering 

Suppose now that one is interested in scattering from a single dielectric sphere S 
of radius R with permittivity ei, permeability and Hclmholtz parameter k\ = 
•\/a; 2 ei/ii, in response to an incoming field (E m ,H m ). The external medium 
is assume to have permittivity e , permeability fi , and Helmholtz parameter 
fco = ^/w 2 eo/io Then the scattered field can be represented by Q with k = ki 
in (|6| for (r, 6, 4>) inside S and by Q with fc = fco in ^ for (r, 0) outside S 1 . 

Let us denote by Eo,Ho the scattered field in the exterior domain and by 
Ei, Hi the scattered field inside S. Then 

E (x) = ««,™ V X V X ( x </>n°m) + iuJ ^0 E & «>™ V X ( X ^m) 

H (x) = Y V™ V X V X ( X ^n%) " iw£ E a «^ iV X ( X ^™) 

where <^ m (x) = </>, fe i:m [r, (9, 0] = /i„(fcr)r„ m (6», 0) and 

Ei(x) = ^c„, m V x V x (x^ m ) + iw/tx^dn^V x (xi/>^ m ) 

n,m n,m 

H i( x ) = E d ">™ V x V x ( x V^m) - iuti ^ C ™,™ V X ( X ^m) 

where V* m (x) = < ro [r, 0, 0] = j n (fcr)F™(0, 0). 

We may also expand (E m ,H m ) in terms of spherical harmonics on the sur- 
face of S: 

E m (x) = a n,™^ X V X (x^ m ) + UJfio E ^n,™V X (x^° m ) 

n.m n,m 

H m (x) = £ ^ n , m V x V x (xV^J - iu;e £ a «- V x ( x ^n?m) 

n.m n.m 

All of the spherical harmonic modes uncouple for fixed n, m, allowing for the 
determination of (a n , m , b n ^ m , c„, m , d„, m ) from the data (a„, m , /3 n ,m) by applying 
the interface conditions ([3]). After some algebra (see, for example, [3J [H]), one 
obtains two uncoupled linear systems of the form 



H n (k R) -J n {hR) \ / \ _ ( -J n [k R)a n 

.m 

(k R) -exjnihR) J \ 

H n (k R) -J n {k x R) \ ( b n , m \ _ f -J7i(koR)l3n,m 
(k R) -nxjnikxR) J \ d 

n . r 



(7) 
(8) 
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where H n (z) = [h n (z) + zh' n (z)}, J n (z) = [j n {z) + zj' n (z)}. 

Definition 2.1. The mapping from incoming coefficients (a nim , ft n ,m) to the 
outgoing coefficients (a„ !m ,6„. m ) is referred to as the scattering matrix and 
denoted by S. 

2.3 Perfect conductors 

If the sphere S is a perfect conductor, the corresponding boundary conditions 
are that the tangential components of the total electric field are zero [T71 [53] : 

n x E tot = n x E = -n x E m . (9) 

In that case, the interior field is identically zero and the scattered matrix is 
given by 

a n ,m = -(Jn(k R) / H n (k R))a n , m 

b n ,m = -{jn{k Q R)/h n (k R))l3 ntm (10) 

3 Scattering from multiple spheres 

Suppose now that one is interested in scattering from M disjoint dielectric 
spheres, where each sphere Si has radius Ri and fc; = \Juj' 2 £ihi. The external 
medium and incoming field are as above. Then, the incoming field can be 
represented on the surface of Si by the expansion 

Ef = E<™ V x V x ^nU^i)) + ^Mo5>U V x ( x 'V^ m (x/)) 

H7 = X V X (*rtnU*l)) - ^o]>>U V X ( X ^%( X 0), 

n.m n.m 

while the scattered field in the interior of Si can be represented by the expansion 
E; = E c «, m V x V x (x^ m (x/)) + iw W d «,m V x ( x ^! m (x/)) (11) 

n,m n.m 

Hi = d n,mV X V X ( x iV4 ! , m (x0) - Yl C U V X ( x ^n,m( x *))- (12) 

n,m n,m 

Here, ip^mix-i) = jn(kri)Y™(9i,(j)i) is computed in terms of the spherical coor- 
dinates (ri, 0i,4>i) of a point x; with respect to the center of Si. 

The scattered field in the exterior of all the spheres can be represented by a 
sum of outgoing expansions, one centered on each sphere. 

M M 



Eo = EE ffl U V x V x (x^° ro (xO) + ^o^^C V x (x,^° m (x z )) 

l—\ n,m 1=1 n,m 

M M 

Ho = X!XX,m V x V x ( Xi ^ m (x0) - zweoE^aLmV x (x ; <^ m (x,)). 



Z=l n,m l — l n,m 
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For a point x exterior to all spheres, the function <^° m (x;) = h n (k ri)Y^ a (9i,4>i), 
where x; = {ri,6i,<j>i), the latter being the spherical coordinates of x with respect 
to the center of Si. The coefficients {a l n , m ,b l nml c l nm ,d l nm ) are all unknowns. 
They are determined by a linear system that imposes the dielectric interface 
condition ^ on each sphere boundary. Unlike the case of a single sphere, 
however, it is no longer trivial to solve for these unknowns, since the incoming 
field experienced on each sphere is due, not only to the known incoming field 
(E m ,H m ), but to the field scattered by all the other spheres. This results in a 
dense linear system involving all of the unknowns, whose solution accounts for 
all of these multiple scattering interactions. 

3.1 Translation operators for multiple scattering 

Fortunately, the outgoing Debye expansion on sphere Sj can be analytically 
converted to an incoming expansion on sphere Si for I =/= j. 

Lemma 1. Let the outgoing expansion from sphere Sj be given by 

E = <™ V X V X ( X Atn( X j)) + ^0 6 »,™ V X ( X ^n°m( X j)) 

H = Yl X V X ( X i^m( X i)) - iuJ6 ° J2 <™ V X ( x ^nV( X j))- 

n,m n,m 

Then, the corresponding field induced on the surface of sphere Si can be repre- 
sented in the form 

K = Yl T&» V X V X ( X ^n?™( X 0) + £ ^' m V X (X^£> m (x0) 

H[, = £ <™ V x V x ( x ^%(*0) ~ »wco E x (x/^? m (x/)). 

n,m n,m 

We denote the mappings from the {a? n m \ and {b 3 n m } coefficients to the {j^ m } 
and {S^ m } coefficients by T°j 7 , Tj'i , T?f , andTh , respectively. Each of these 
mappings depends on the vector from the center of sphere Sj to sphere Si and 
the parameters (jj,Q,eo,oj). 

For convenience, we will sometimes denote vectors of coefficients such as 
{ a n,m} by a- 7 - The individual components of a translated vector such as T?j <P 
will be denoted by [T?f '<P] n>m . 

Remark 3.1. The formulae for the translation operators T°J 7 , T^'j 7 , Tjf, and 
Tjf are rather involved JSj 1 1 6 A \21f . If the expansions are truncated at n = p 
terms, there are 2(p+ f) 2 nonzero coefficients in both the outgoing (a J n m ,V n m ) 
and incoming ("f^ l m , <^, m ) representations. Each translation operator is dense 
and, therefore requires 0{p A ) operations to apply. More efficient schemes f$l\lH$ 
reduces the cost to 0(p 3 ), while the diagonal- form of the EMM £3 \2^j reduces 
the cost to 0{p 2 \ogp) for well-separated spheres in the high-frequency regime. 
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Let us now assume that all outgoing and incoming expansion are truncated 
at n = p terms. The choice of p is determined by accuracy considerations. It 
must be sufficiently large to resolve the E and H fields on each sphere surface 
to the desired precision. 

Using the preceding lemma, the total field immediately exterior to sphere Si 
can be written in the form 

M 

E l = E« + Y. I7"' ! + rjy^k™ V x V x (x^ m (x ( )) 
J=l 

M 

+ iwoY^lTffa) + T b j9] n , m V x (x^ m (xO) (13) 



3=1 
3*1 



V 



( x 0) + * w Mo^ 

V x (x^ 

( x 0) 

M 

H< = Hf + £[7^" + V x V x (x,^(x,)) 

i=i 
A/ 

- iu,e $><y<z3 + ijf V x (x^ ro (*0) (14) 

j=i 

P P 

+ XX™ V X V X ( X ^nV( X 0) " *Weo5^°n.w V X ( X ^n%( x 0) ■ 



The first terms in the preceding expressions for Eq, Hq account for the incoming 
field, while the next two terms account for the scattered field coming from all 
other spheres. The last two terms in each expression account for the fields being 
scattered by Si itself. 

It is now clear how to apply the interface conditions (|3|). We simply equate 



the tangential components of Eq,Hq defined in ( 13 ) , (I14J) with the tangen- 



tial components of the interior representations (E/,H;) defined in (11 1,(12). 
This yields a dense linear system of dimension 4M(p + l) 2 for the coefficients 
(a l n m , b l n m , c l n m , d l n m ) . We will refer to this system as the multiple scattering 
equations. Writing the equations out explicitly is not especially informative, 
and we omit it. 



Remark 3.2. The scattering matrix S (Definition 2.1) allows for the elimina- 
tion of the interior variables (c n ml d l n m ), so that the one can solve a modified 
system of dimension 2M(p + l) 2 for the coefficients (a l n m , b l n m ) describing the 
exterior field alone. 
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M 



a 




) 



S 



(15) 



A I 



n,m 



\ 



3 = 1 



J 



It is worth emphasizing that the multiple scattering equations are hardly 
new. There is a vast literature on the subject, which we do not seek to review 
here. We refer the reader to the textbooks [2J |3J [T3J [TS1 [3T] and the papers 



3.2 Iterative solution of the multiple scattering problem 
for a system of spheres 

We will solve the multiple scattering equations iteratively, using GMRES [2"5] 
with a block diagonal preconditioner, each block corresponding to the unknowns 
on a single sphere. In applying the preconditioner, we simply invert each of the 
M diagonal blocks, which corresponds to solving the single sphere scattering 
problem described in section |2.2| Since all M spheres interact, however, the 
system matrix is dense. Each matrix-vector multiply in the iterative solution 
process, if carried out naively, would require 0(M p ) work. 

In order to accelerate the solution procedure, the wideband fast multipolc 
method (FMM) [5] can easily be modified to reduce the cost to 0(Mp 3 ) work 
per iteration. This is discussed in the context of acoustic scattering in [T31 H] • 
Since the literature on FMMs is substantial, we omit a detailed discussion of 
the technique, but present results in section [6] 

4 Scattering from an arbitrary inclusion 

Suppose now that instead of a sphere, we are given a smooth inclusion (or set of 
inclusions) D\ with permittivity e% and permeability p,\ embedded in the same 
infinite medium as above. We will suppose further that D\ can be enclosed 
in a sphere S\ (Fig. [2]). As before, at the material interface, the conditions 
to be satisfied are ((3j). The Debye-Lor enz-Mie formalism cannot be applied in 
this case, and attempts to do so (called the T-matrix method) suffer from ill- 
conditioning when Di is sufficiently non-spherical. We, therefore, turn to the 
standard representation of electromagnetic fields in general geometries, based 
on the vector and scalar potentials and anti-potentials [551 ■ 
The vector potential in domain I (1 = 0, 1) is defined by 



men]. 
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Figure 2: A pair of triangulated ellipsoids define a bounded domain D\ that lies 
with an enclosing sphere S\. The scattering matrix for D\ will be created on S\ 
and used to represent the exterior field. 



and gi(x) — e jfei " x "/|| x ll with k\ — ^Jw 2 ej^j. When the argument of the square 
root is complex, fc; is taken to lie in the upper half-plane. We define the vector 
anti-potential in domain I by 



Aj(x)=ej/ gi{x-y)K[ds. 



From these, we may write 



where 



1 



= + ICjAf V X A; 

H; = — V X A; - Vipi + icoAi. 



tujeifii 
1 

iujeifii 



-V- A, 



-V- A, 



As written above, we have twelve degrees of freedom at each point P € 
dDi, namely the three Cartesian components of J , Ji,K ,K l5 but only four 
boundary conditions (the continuity of the tangential components of E and H) . 
We will assume, however, that the functions Jo, Ji, Ko, Ki are surface currents 
and that the following linear relations hold 



£o 



Ji 



Kn 



/'i 



Ki. 



This leaves four degrees of freedom. Imposing the conditions ([3| on Ji,Ki 
results in Muller's integral equation [22] , a resonance-free Fredholm equation of 
the second kind. 
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Figure 3: In the simplest geometric model, the surface of the scatterer dDi is 
approximated by a collection of flat triangles, defined by the locations of its three 
vertices in M 3 . On each triangle, there are two two linearly independent tangent 
directions ti and t 2 . The unknown electric and magnetic currents J and K on each 
triangle are defined by Jiti+,72t2 and fciti+fc2t2, respectively, and the electromag- 
netic fields are evaluated at the triangle centroids. For higher order accuracy, each 
quadratic surface patch is specified by six nodes: the three triangle vertices and 
three additional points, one on each curved triangle side. Three "support nodes" 
x^x^x 3 are then selected in the interior of each patch. Our representation for 
J and K at each support node x J is of the form j\t\ + and k\t\ + k\t\, 
where t^,t 2 are linearly independent tangent vectors at x\ The support nodes are 
also the points where we evaluate the electromagnetic fields and impose interface 
conditions. 

In more detail, using the facts that 

V x x (0,(x - y) K(y)) = V x9l x K(y) , 
ax b x c = b(a ■ c) — c(a ■ b) , 

and, for y € dD\, 

A™ / 7^f"( x ~ y ) <j (y) ds y = l a (yo) + / 7?r L (yo - y) ^(y)d» y 

x€ n° JdD 1 on ya I J dDl 0n yo 

/ 7^ L (*-y) a (y) ds y = -l a (y°)+ ( f Sr-(yo-y)v(y)ds y , 

xe™ JdDt on ya I J 9Dl On ya 

we obtain the following coupled set of equations. 



n x W 



— [eoMo.9o - eiA*i5i] (n x Ji)ds y 

£ 1 JdD t 

H — — n x / [VV.9o - VV51] Ji ds y 
Mo/ f^-^)(n.K0^ 



(16) 



+ 



1 



1 



2e 2ei 



Ki + Hq 



1 dg 



1 dg 



dDi \ n i^a dn n ei dn 



- I ds y 
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n x H' m = 



- [Moeoffo - eiMiffi] (n X K x ) ds y 

1 JdDi 

-n x / [VV.go - VV.gi] K x ds y 

1 JdDt 



I 

UJjJL 



(17) 



-eo/ p°_YM (n . Ji} 



2/i 2fiiJ 1 °7ar> 1 V e iMo^ e [ii dn J 1 y ' 

Because the Miiller equation is a second kind Fredholm equation, the order 
of accuracy of the solution is that of the underlying quadrature rule. For first 
order accuracy, we assume Ji and Ki are piecewise constant current densities 
on a flat triangulated surface. For second order accuracy, we assume Ji and Ki 
are piecewise linear current densities on a piecewise quadratic surface with each 
curved triangle defined by six points (Fig. [3|. 

For each discretization node, we evaluate the relevant electromagnetic field 
component using a mixture of analytic and numerical quadratures on each tri- 
angle. More precisely, we use the method of singularity subtraction - computing 
integrals analytically for the kernel 1/r and its derivatives and using numerical 
quadrature for the difference kernel \e %kr — l]/r, which is smoother. This results 
in a complex linear system of dimension AN x AN for first order accuracy and of 
dimension 127V x 127V for second order accuracy, where N denotes the number 
of triangles. For small N, say N < 1000, one can use direct LU-factorization to 
solve the linear system. For larger values of N, iterative solution with FMM- 
accelcration becomes much more practical [5J 15] . 

4.1 The scattering matrix for Di 

Suppose now that we are interested in scattering from the two ellipsoids Di 
shown in Fig. [2] due to an incoming field which is regular in the enclosing 
sphere Si- Such an incoming field can be expanded within S% in the form 

E m (x) - a 'M»V x V x (x^° m ) + iw/*, x ( x ^»%) 

H"( X ) - / 3 -" V X V X ^n°m) - iUtO «™^ V X i^n°m) , 



as in Section 2.2 Each (vector) spherical harmonic modes, corresponding to a 
single a„ jm or /3 njm , defines a particular incoming field on D\. More precisely, 
we can solve the Miiller equation for a right-hand side obtained by setting the 
incoming field to be 

Ei"(x) = V x V x (x^° ), Hi" (x) = - lW e V x (x^° ) , (18) 
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corresponding to setting a fixed a n ^ m = 1 and all other coefficients to zero. 
Similarly, we can set the incoming field to be 

E|%, ro (x) - + ^MoV x (x^ m ), H» >ro (x) = V x V x (x^° m ) (19) 

corresponding to setting a fixed /3 n , m = 1 and all other coefficients to zero. 
We can then store either the electric and magnetic currents j\' n ' m , Kj' n,m or 
j2, n,77^ j^2,n,m m( j ucec i by these (unit) incoming fields or just convert these 
currents to the coefficients of the outgoing (scattered) fields: 

Ef„, m (x) = J2 * V X ^n,,m>) + ™W> £ X (x^,) 

n' \ra' n' ,m' 

H?:„, m (x) = £ 6##V x V x (x<^ m ,) - iW£Q £ o$£V x (x<^ m ,) 
and 

Ei:„, m (x) = £ <#£?V x V x (x^ >m/ ) + iww, £ ft^V X (x^* m ,) 

H^, m (x) = £ 6^V x V x (x^ m ,) - iW£o £ <#£?V x (x^ m ,). 

The formula for converting the currents jj :n ' m , K}' n,m to the coefficients can be 
obtained by orthogonal projection of the induced field on the enclosing sphere 

By superposition, an incoming field defined by the vector of incoming coef- 
ficients {a n ,m, /?n, m } results in a scattered field of the form 

E sc (x) = £ «n',m'V X V X (x^ m ,) + lu ^ Q &*',m'V X (x^? >m ,) 

n',m ; n\m' 

H sc (x) = £ ftn'.tn'V X V X (X0*? im ,) - iweo £ «™',™' V X ( x ^n',m') ' 
n' ,m' n' ,m' 

with the coefficients of the scattered field given by 

El,n,m i /O 2.n,m 

n,m 

Et l,n,m l' 



n,m 



2,n,m 
n' ,m' " 



The matrix mapping the incoming to the scattered coefficients is referred to as 
the scattering matrix for the structure Di- 

Fixing the order of the expansions above at p, there are 4p 2 possible basis 
functions that span the space of all possible incoming fields. We must, there- 
fore solve 4p 2 Miiller integral equations on the detailed geometry defining D\. 
To store the currents induced by each incoming field requires 0(Np 2 ) memory, 
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where N denotes the number of degrees of freedom used in the discretization 
of the integral equation. The scattering matrix itself requires storing 0(16p 4 ) 
complex numbers. While somewhat expensive, this is a pre-computation step, 
in anticipation of simulating microstructures with thousands or millions of inclu- 
sions of the same identical shape, but well enough separated that the scattering 
matrices are accurate. 

5 Multiple scattering from well-separated non- 
spherical inclusions 

Once the scattering matrix is known, the solution to the full Maxwell equations 
for geometries with N inclusions (N — 200 in Fig. [lj can be turned into 
a multiple-scattering problem based only on the enclosing spheres. That is, 
the inclusions can be replaced by their scattering matrices and the multiple 
scattering method of section [3] can be used with trivial modifications. 

There are two distinct advantages to be gained here. First, we have re- 
duced the number of degrees of freedom from, say, 5,000 or 10,000 unknown 
current density values per inclusion to, say, 400 expansion coefficients. Just 
as important, however, is that we have precomputed the solution operator for 
each inclusion in isolation, so that the linear system we solve by iteration on 
the multi-sphere system is much more well-conditioned. Further, the FMM re- 
duces the cost of each iteration from 0(N 2 ) to 0(N log N) and is particularly 
efficient here, since the complicated quadratures on triangulated surfaces have 
been subsumed into the precomputation step. 

The principal limitations of the method are 1) that some modest separation 
distance between inclusions is required and 2) that the bookkeeping becomes 
a bit awkward if more than a few distinct nanoparticle types are allowed. In 
many experimental settings, both conditions are satisfied. 

It is worth noting that the method of this paper can be viewed as a reduced 
order model for the scattering problem. In broad terms, the idea is not new and 
there is substantial activity in this area in both electromagnetics and other fields 
(see, for example, [3]). It is also worth noting that the method is "rigorous" 
in the sense that the error is determined in a straightforward manner by the 
accuracy of the Miiller integral equation solver and the order of expansion of 
the scattering matrix. It fails (or needs local modification) if and only if two 
enclosing spheres intersect. 

6 Numerical Examples 

As discussed in section |4j the Miiller integral equation is an effective method for 
determining the scattering matrix from a dielectric inclusion of arbitrary shape. 
To illustrate its performance, we consider the geometry in Fig. [2j consisting of 
a pair of ellipsoids triangulated with piecewise quadratic triangles on which we 
allow piecewise linear current densities. Each triangle has three nodes with two 
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degrees of freedom for each current (electric and magnetic) at each point, result- 
ing in a complex linear system of dimension 2160 x 2160. (All calculations and 
timings reported in this section have been carried out using a 12-core 2.93GHz 
Intel Xeon workstation.) LU factorization requires 3.5 seconds, and the subse- 
quent solution requires 0.1 seconds for each possible incoming mode. With 720 
triangles, the linear system has dimension 8640 and with 2880 triangles, there 
are 34,560 degrees of freedom. These require 106 and 2,620 seconds to factor, re- 
spectively. The solution times for each incoming mode are 3.52 and 87 seconds, 
respectively. We could accelerate these solution times using fast multipole-based 
codes (or any of a variety of other "fast" algorithms) , but we view this cost as 
an initialization step and the CPU times are acceptable. The errors are of the 
order 10 -3 , 10 -4 , and 10 -5 for the successively finer discretizations, somewhat 
better than the expected second order convergence. 

To illustrate the performance of the FMPS algorithm, we consider a 21 x 
21x2 array of scatterers, each consisting of an ellipsoid pair with a scattering 
matrix derived from the Miiller integral equation of order p = 3. Using the 
same 12-core 2.93GHz Intel Xeon workstation, the time required was about 2 
seconds per iteration, with six iterations required for GMRES to converge to 3 
digits. The "slow" multiple scattering (SMPS) approach, without fast multipole 
acceleration, required about 7 seconds per iteration. For a 21 x 21 x 4 array, 
the cost was about 6 seconds per iteration (28 seconds for SMPS) and for a 
21 x 21 x 8 array, the cost was about 23 seconds per iteration (108 seconds for 
SMPS). For a 21 x 21 x 16 array (14,112 ellipsoid pairs), the cost was about 59 
seconds per iteration (440 seconds for SMPS). 

The reason for the modest speedup of the FMPS over the SMPS approach is 
that the number of spheres is still rather small. For one million scatterers, the 
speedup factor would be about 1000. Careful readers may note that the FMPS 
scaling appears worse than 0(N log N) in successively doubling the simulation 
from a 21 x 21 x 4 array to a 21 x 21 x 8 array to a 21 x 21 x 16 array. For those 
familiar with the FMM, the short explanation is that the "near neighbor" cost 
is not yet in the asymptotic regime in the first two cases. Timings extrapolated 
from the last case are accurate for any volume-filling distribution. 

Finally, we illustrate the use of the FMPS algorithm in carrying out fre- 
quency scans for (a) one ellipsoid pair with the long axis oriented parallel to the 
(linearly polarized) incoming electric field, (b)one ellipsoid pair with the long 
axis oriented parallel to the (linearly polarized) incoming magnetic field, or (c) 
four pairs of randomly oriented ellipsoid pairs (Fig. [4|. 

7 Conclusions 

The method introduced in this paper (fast multi-particle scattering) combines a 
highly accurate integral equation solver with multiple scattering theory, in order 
to permit the solution to the full Maxwell equations in configurations typical of 
engineered composites (metamaterials). We assume that the geometry consists 
of a large number of inclusions embedded in a homogeneous background. While 
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Figure 4: The top row shows a frequency scan of the real and imaginary parts 
of electric polarization vector (left), the real and imaginary parts of magnetic po- 
larization vector (middle), and the scattering (right, upper curve) and absorption 
(right, lower curve) for one ellipsoid pair with the long axis oriented parallel to 
the (linearly polarized) incoming electric field. The second row shows the same 
computed quantities for one ellipsoid pair with the long axis oriented parallel to 
the (linearly polarized) incoming magnetic field. The third row shows the same 
computed quantities for four pairs of randomly oriented ellipsoid pairs. 
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we have only included a single type of inclusion geometry in our examples above, 
it is clear that the method can easily be applied to permit several such types, so 
long as there is a modest separation between inclusions. FMPS is enormously 
faster than a full FMM-based solver using the full discretization of the geometry. 
With 14,112 ellipsoid pairs (the largest example in the preceding section), this 
would require about 30 million degrees of freedom, many minutes per iteration, 
and many more iterations. 

In its present form, the method cannot be used for tightly packed configu- 
rations, which will require more elaborate compression schemes [llj . It does, 
however, permit workstation-based simulation with millions of inclusions. We 
are currently working on extending the method so that it can handle inclusions 
embedded in a layered medium. 
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